%% 设置随机种子以获得可重复的结果
clc,clear;
close all;
warning off;
timetable=readtable("time_series_data.csv");
time=timetable(:,1);
value=timetable(:,8);
value=table2cell(value);
value=cell2mat(value);
% 分割训练集和测试集
trainSize = round(0.8 * length(value));
trainData = value(1:trainSize);
testData = value(trainSize+1:end);
%% 数据图
figure(1);
plot(value);
title('Time Series Data');
%% 季节性分析
% 绘制周期图
ts=timeseries(value);
figure(2);
periodogram(ts.Data);
title('周期图');

% 频谱分析
[f, Pxx] = periodogram(ts.Data, [], [], 1); % 假设采样频率为1(每天)
[maxP, idx] = max(Pxx); % 找到功率谱密度的最大值
seasonality = round(1 / f(idx)); % 计算对应的周期
s=seasonality;
% 输出季节性周期
disp(['推测的季节性周期为: ', num2str(seasonality)]);
%% d值的确定（非季节性差分数）
d=0;
flag=adftest(trainData);%平稳性检验
while flag==0 
    d=d+1;
    hk=diff(trainData,d); 
    flag=adftest(hk);
end
%% BIC 准则定阶 AIC 准则定阶(用AIC的线比较好BIC简单)
maxLags = 5;
BIC = NaN(maxLags+1,maxLags+1);
AIC = NaN(maxLags+1,maxLags+1);
LOGL=zeros(maxLags+1,maxLags+1);
PQ=zeros(maxLags+1,maxLags+1);
n=length(trainData);
for p = 0:maxLags
    for q = 0:maxLags
        try
            % 定义ARIMA模型
            mdl = arima(p, d, q);
            
            % 使用evalc捕获并抑制输出
            [l,s,EstParamCov,LOGL] = evalc('estimate(mdl, trainData, "Display", "off");');
            
            % 提取估计参数之间对数似然函数最优值
            LogL = LOGL;
            
            % 计算有多少个估计参数
            numParams = sum(any(EstParamCov));
            % 计算AIC和BIC
            [AIC(p+1, q+1), BIC(p+1, q+1)] = aicbic(LogL, numParams, n);
        catch
            % 如果估计失败，则继续下一个组合
            continue;
        end
    end
end
N = min(min(BIC)); [p1,q1] = find(N==BIC);
p1=p1-1;q1=q1-1;
% 画热度图来表示 BIC 数值的分布
figure(3)
heatmap(BIC/100);
xlabel("AR")
ylabel("MA")
title('Optimal ARIMA(BIC) ')
M = min(min(AIC)); [p2,q2] = find(M==AIC);
p2=p2-1;q2=q2-1;
% 画热度图来表示 AIC 数值的分布
figure(4)
heatmap(AIC/100);
xlabel("AR")
ylabel("MA")
title('Optimal ARIMA(AIC) ')

%% ARIMA模型
Mdl =  arima(p2,d,q2);
EstMdl = estimate(Mdl, trainData, "Display", "off");
numSteps = length(testData);%检验
[YF, YMSE] = forecast(EstMdl, numSteps, 'Y0', trainData);
num = 10;%预测数
[YF1, YMSE1] = forecast(EstMdl, num, 'Y0', value);
nk=numel(value);
%预测
disp('Forecasted Values:');
disp(YF1);
disp('Forecast Mean Squared Error:');
disp(YMSE1);

lower = YF -1.96*sqrt(YMSE); %95置信区间下限
upper = YF + 1.96*sqrt(YMSE); %95置信区间上限
lower1 = YF1 -1.96*sqrt(YMSE1); %95置信区间下限
upper1 = YF1 + 1.96*sqrt(YMSE1); %95置信区间上限

%% 误差与预测
forecastError = testData - YF;
ME=sum(forecastError)/numSteps;  %平均误差
REME=sqrt(sum(forecastError.^2)/numSteps);   %均方相对误差
Emax=max(forecastError);   %最大误差
figure;
hold on;
plot(1:length(trainData),trainData, 'b');
plot(length(trainData):length(value),[trainData(end);testData], 'r');
plot(length(trainData):length(value),[trainData(end);YF], 'k');
plot(length(trainData):length(value),[trainData(end);lower],'r:','LineWidth',2);
plot(length(trainData):length(value),[trainData(end);upper],'r:','LineWidth',2);
xlabel('Time');
ylabel('Value');
legend('Forecasted Values','95% 置信区间');
hold off;

figure;
hold on;
h0=plot(value, 'b');  % 原始数据
h2=plot(numel(value):numel(value)+num,[value(end);lower1],'r:','LineWidth',2);
plot(numel(value):numel(value)+num,[value(end);upper1],'r:','LineWidth',2);
h1=plot(numel(value)+(1:num), YF1, 'k');  % 预测数据
line([nk,nk+1],[value(nk),YF1(1)]);%修线
xlabel('Time');
ylabel('Value');
legend([h0,h1,h2], 'Forecasted Values','95% 置信区间');
hold off;

%% 检查数据的特性(ACF/PACF(P,Q)，残差)

% 计算残差
residuals = infer(EstMdl,trainData);
% 绘制图
hold on;
figure(6);
subplot(3,1,1);
plot(residuals);
title('Residuals');
xlabel('Time');
ylabel('Residual');
subplot(3,1,2);
autocorr(value);
title('ACF');
subplot(3,1,3);
parcorr(value);
title('PACF');

%% 检验

% Q-Q图
figure;
qqplot(residuals);
title('Q-Q Plot of Residuals');
% Ljung-Box检验(h=1,存在自相关性)
[h, pValue, stat, cValue] = lbqtest(trainData, 'Lags', 10);
disp('Ljung-Box Test:');
disp(['h = ', num2str(h)]);
disp(['pValue = ', num2str(pValue)]);
disp(['stat = ', num2str(stat)]);
disp(['cValue = ', num2str(cValue)]);
% 解释结果
if h == 0
    fprintf('不能拒绝零假设，时间序列没有显著的自相关性（白噪声）。\n');
else
    fprintf('拒绝零假设，时间序列存在显著的自相关性。\n');
end
